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Abstract: Cumulative link models have been widely used for ordered categorical 
responses. Uniform allocation of experimental units is commonly used in practice, 
but often suffers from a lack of efficiency. We consider D-optimal designs with 
ordered categorical responses and cumulative link models. For a predetermined set 
of design points, we derive the necessary and sufficient conditions for an allocation 
to be locally D-optimal and develop efficient algorithms for obtaining approximate 
and exact designs. We prove that the number of support points in a minimally 
supported design only depends on the number of predictors, which can be much 
less than the number of parameters in the model. We show that a D-optimal 
minimally supported allocation in this case is usually not uniform on its support 
points. In addition, we provide EW D-optimal designs as a highly efficient surrogate 
to Bayesian D-optimal designs. Both of them can be much more robust than 
uniform designs. 

Key words and phrases: Approximate design, exact design, multinomial response, 
cumulative link model, minimally supported design, ordinal data. 

1. Introduction 

In this paper we determine optimal and efficient designs for factorial 
experiments with qualitative factors and ordered categorical responses, or 
simply ordinal data. Design of experiment with multinomial response, 
and ordered categories in particular, is becoming increasingly popular in 
a rich variety of scientific disciplines especially when human evaluations 
are involved (Christensen, 2015). Examples include wine bitterness study 
(Randall, 1989), potato pathogen experiments (Omer et al., 2000), radish 
seedling’s damping-off study (Krause et ah, 2001), polysilicon deposition 
study (Wu, 2008), beef cattle research (Osterstock et ah, 2010), and toxic¬ 
ity study (Agresti, 2013). 

This research is motivated by an odor removal study conducted by the 
textile engineers at the University of Georgia. The scientists studied the 
manufacture of bio-plastics containing odorous volatiles, which need to be 
removed before commercialization. For that purpose, a 2 2 factorial exper¬ 
iment was conducted using algae and synthetic plastic resin blends. The 
two factors are types of algae (aq: raffinated or solvent extracted algae 
(—), catfish pond algae (+)) and synthetic resins (x 2 : polyethylene (—), 
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polypropylene (+))• The response Y has three ordered categories: serious 
odor (j = 1), medium odor (j = 2), and almost no odor (j = 3). Follow¬ 
ing the traditional factorial design theory, a pilot study with equal number 
(10 in this case) of replicates at each experimental setting was conducted, 
known as a uniform design. The results are summarized in Table 1, where 
Dij represents the number of responses falling into the jth category under 
the ith experimental setting. As demonstrated later (Section 4), the best 
design identified by our research could improve the efficiency by 25% with 
only three experimental settings involved. 


Table 1: Pilot Study of Odor Removal Study 


Experimental 

Factor 

level 

Summarized responses (Y, odor) 

setting 

Algae 

Resin 

Serious 

Medium 

No odor 

i 

Xi 

X2 

yn 

Vii 

ya 

1 

+ 

+ 

2 

6 

2 

2 

+ 

- 

7 

2 

1 

3 

— 

+ 

0 

0 

10 

4 

- 

- 

0 

2 

8 


For such kind of ordinal response Y with J categories and d predictors 
x = (xi,..., Xd) T , the most popular model in practice was the proportional 
odds model (also known as cumidative logit model , see Liu and Agresti (2005) 
for a detailed review). McCullagh (1980) extended it into the cumulative 
link model (also known as ordinal regression model) 

g (P(Y < j | x)) = 0j =- /3 T x, J = (1.1) 

where g is a general link function. The proportional odds model is a special 
case when g is the logit link. Nevertheless, link functions other than logit 
may fit some experiments better. Examples include the complementary 
log-log link for the polysilicon deposition study (see Example 4.6) and the 
cauchit link for the toxicity study (see Example 5.9). In this paper, we 
adopt the cumulative link model (1.1). 

When there are only two categories (J = 2), the cumulative link model 
(1.1) is essentially a generalized linear model for binary data (McCullagh 
and Nelder, 1989; Dobson and Barnett, 2008). For optimal designs un¬ 
der generalized linear models, there is a growing body of literature (see 
Khuri et al. (2006), Atkinson et al. (2007), Stufken and Yang (2012) and 
references therein). In this case, it is known that the minimum number 
of experimental settings required by a nondegenerate Fisher information 
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matrix is d + 1, which equals the number of parameters (Fedorov, 1972; 
Yang and Mandal, 2015). A design with the least number of experimental 
settings, known as a minimally supported design , is of practical significance 
with a specified regression model due to the cost of changing settings. It is 
also known that the experimental units should be uniformly assigned when 
a minimally supported design is adopted for binary response or under a 
univariate generalized linear model (Yang and Mandal, 2015). 

When J > 3, the cumulative link model should rather be treated as a 
special case of the multivariate generalized linear model (McCullagh, 1980). 
The relevant results in the optimal design literature are meagre and re¬ 
stricted to the logit link function (Zocchi and Atkinson, 1999; Perevozskaya 
et ah, 2003). In this paper, we obtain both theoretical results and efficient 
algorithms for general link functions and reveal that the optimal designs 
with J > 3 are quite different from the cases with J = 2. We prove that 
the minimum number of experimental settings is still d+1 but strictly less 
than the number of parameters d + J — 1 (Theorems 2.3 and 2.4). This 
counter-intuitive result is due to the multinomial-type responses: from a 
single experimental setup, the summarized responses have J — 1 degrees 
of freedom, requiring fewer number of distinct experimental settings in a 
minimally supported design. Because of the same reason, the allocation 
of replicates in a minimally supported design is usually not uniform either 
(Section 5), which is different from the traditional factorial design theory. 

Similar to the cases of generalized linear models, the information matrix 
under cumulative link models depends on unknown parameters. In the lit¬ 
erature, different approaches have been proposed to solve the dependence of 
optimal designs on unknown parameters, including local optimality (Cher- 
noff, 1953), Bayesian approach (Chaloner and Verdinelli, 1995), maximin 
approach (Pronzato and Walter, 1988; Imhof, 2001), and sequential proce¬ 
dure (Ford et ah, 1989). As pointed out by Ford, Torsney, and Wu (1992), 
locally optimal designs are not only important when good initial parame¬ 
ters are available from previous experiments, but can also be a benchmark 
for designs chosen to satisfy experimental constraints. In this paper, we 
mainly focus on locally optimal designs. For situations where local val¬ 
ues of the parameters are difficult to obtain but the experimenter has an 
idea of the range of parameters with or without a prior distribution, we 
recommend EW optimal designs, where the Fisher information matrix is 
replaced by its expected values (Atkinson et ah, 2007; Yang, Mandal, and 
Majumdar, 2016). We compare Bayesian D-optimal designs (Chaloner and 
Verdinelli, 1995) with EW D-optimal designs for ordinal data. As a sur¬ 
rogate for Bayesian designs, EW design is much easier to find and retains 
high efficiency with respect to Bayesian criterion (Section 6). 
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Among various optimal design criteria, D-optimality, which maximizes 
the determinant of Fisher information matrix, is the most frequently used 
one (Zocchi and Atkinson, 1999) and often performs well according to other 
criteria (Atkinson et ah, 2007). In this paper, we study D-optimal designs. 

In the design literature, one type of experiments deal with quantitative 
or continuous factors only. Such a design problem includes identification of 
a set of design points {x;}j=and the corresponding weights {pi}i= 

(see, for example, Atkinson et al. (2007) and Stufken and Yang (2012)). For 
this type of optimal design problems, numerical algorithms are typically 
used for cases with two or more factors (see, for example, Woods et al. 
(2006)). Another type of experiments employ qualitative or discrete factors, 
where the set of design points is predetermined and only the 

weights {pi}i =are to be optimized (see, for example, Yang and Mandal 
(2015)). One connection between the two types of designs is that one can 
pick up grid points of the continuous factors and turn the first-type problems 
into the second. Tong, Volkmcr, and Yang (2014, Section 5) also bridged the 
gap between the two types of problems in a way that the results involving 
discrete factors can be applied to the cases with continuous factors as well. 
In this paper, we concentrate on the second type of design problems and 
assume that {xj}j = i r .. jm are given and fixed. 

This paper is organized as follows. In Section 2, we describe the prelim¬ 
inary setup and obtain the Fisher information matrix for cumulative link 
model with a general link, which generalizes Perevozskaya et al. (2003)’s re¬ 
sult for logit link. We also identify a necessary and sufficient condition for 
the Fisher information matrix to be positive definite, which determines the 
minimum number of experimental settings required. In Sections 3 and 4, 
we provide both theoretical results and numerical algorithms for searching 
locally D-optimal approximate or exact designs. In Section 5, we identify 
analytic D-optimal designs for special cases to illustrate that a D-optimal 
minimally supported design is usually not uniform on its support points. 
In Section 6, we illustrate by examples that the EW D-optimal design can 
be highly efficient with respect to Bayesian D-optimality. We make con¬ 
cluding remarks in Section 7 and relegate additional proofs and results to 
the Supplementary Materials. 

2. Fisher Information Matrix and Its Determinant 

Suppose there are m (m > 2) predetermined experimental settings. For 
the ith experimental setting with corresponding predictors x,; = (xn, ..., 
Xid) T E (d > 1), there are n t experimental units assigned to it. Among 
the rij experimental units, the 7th one generates a response VA which be¬ 
longs to one of J (J > 2) ordered categories. As shown in Example 2.2 
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later, the dimension d of the predictors can be significantly larger than the 
number of factors considered in the experiment, which allows more flexible 
models. 

2.1 General setup 

In many real applications, Vn,..., I/ ni are regarded as i.i.d. discrete ran¬ 
dom variables. Denote 7 = P(I4t = j), where i = 1,..., m; j = 1,-, J; 

and k = 1,..., n* . Let Y Z] = #{k \ Vik = j} be the number of V^’s falling 
into the jth category. Then (Y, u ..., Y tJ ) ~ Multinomial(n*; tth, ..., 7 qj) . 
Throughout this paper, we assume 

Assumption 2.1. 0 < 7 T*j- < 1, i — 1,..., m; j — 1,..., J. 

Denote 7 ^ = P(V ik < j) = irn + ■ ■ ■ + 1 Uj , j = 1,..., J. Based 
on Assumption 2.1, 0 < 7*1 < 7*2 < • • • < 7 «,j-i < 7 ij = 1 for each i = 
1 ,... ,m. Consider independent multinomial observations (l'h, • • •, Yjj),i = 
1 , ..., m with corresponding predictors x 1? ..., x m . Linder a cumulative 
link model or ordinal regression model (McCullagh, 1980; Agresti, 2013; 
Christensen, 2015), there exists a link function g and parameters of interest 
6 *i,..., Oj-i, (3 = {/3i ,..., (3 d ) T , such that 

9(Kij) = 0j - *1 P'y j = Uky ■, J ^:1- (2- 1 ) 

This leads to m( J — 1) equations in d+ J — 1 parameters (/?i ,..., /3 d , 0i,..., 
9j_i). Throughout this paper, we assume 

Assumption 2.2. The link g is differentiable and its derivative g' > 0. 

Assumption 2.2 is satisfied for commonly used link functions includ¬ 
ing logit (log( 7 /(l - 7 )), probit ($ - 1 ( 7 )), log-log (-log(-log( 7 ))), 
complementary log-log (log(—log(l— 7 ))), and cauchit (tan(7r(7—1/2))) 
(McCullagh and Nelder, 1989; Christensen, 2015). Some relevant formulas 
of these link functions are provided in the Supplementary Materials (Sec¬ 
tion S.l). According to Assumption 2 . 2 , g is strictly increasing, and then 
61 < 62 < ■ ■ ■ < Qj- 1 . 

Example 2.1. Consider the logit link < 7 ( 7 ) = log( 7 /(l — 7 )) with two 
predictors and three ordered categories. Model (2.1) consists of 2m equa¬ 
tions g(jij) = 9j — Xnfii — 2 , i — 1,..., m; j — 1,2 and 4 parameters 

(Al, P 2 , 9\, 62 )- Linder Assumptions 2.1 and 2 . 2 , 9\ < 62 ■ □ 

Example 2.2. Suppose the model consists of three covariates X\,X 2 , X:> and 
a few second-order predictors, < 7 ( 7 ^) = dj — xnfdi — x^^ — Xisfds — XiiX^fdu — 
xfiPn — x^ 2 j 322 , where i = 1 ,,m;j = 1,..., J — 1. Then the number of 
predictors is d = 6 . □ 
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Under the cumulative link model (2.1), the log-likelihood function (up 
to a constant) is l(f} u ..., /3 d , 0 1? ..., 6 j_ i) = YZLi Yfj =1 Y H l°g(^j), where 
7 Tij = 7 ij - jij-i with 7 ij = - xf/3) for j = 1,..., J - 1 and 7*0 = 0, 

lij = 1 - 

Perevozskaya et al. (2003) obtained a detailed form of the Fisher infor¬ 
mation matrix for logit link and one predictor. Our theorem below is for 
general link and d predictors with proof relegated into the Supplementary 
Materials (Section S.3). 

Theorem 2.1. Under Assumptions 2.1 and 2.2, the Fisher information 
matrix can be written as 

m 

F = ^n*A* ( 2 - 2 ) 

i= 1 

where A* is a (d + J — 1) x (d + J — 1) matrix 

A{ 1 Aj2 \ _ f (Ci3'is*£it)s=l,...<i;t=l,...,ci ( Xi s Cnfs=l,...,d\t=l,...,J—l 

Aj 2 Aj 3 J y ( Cis^it)s=l,...,J—l‘,t=l,...,d Aj 3 

and Aj 3 is a (J — 1 ) x (J — 1 ) symmetric tri-diagonal matrix with diagonal 
entries Un,..., and off-diagonal entries — 6 * 2 , • • •, — 6 *,j-i when J > 3, 

where e* = - g id _ i ) 2 > 0 with g l} = (. g "^'(^ - xf/3) > 0 for 

j 1, • • •, J 1 and gio g^j 0, Cn 9it\Fit i.9it 9i,t— i) ^"*,t+i(5 , *,t+i 
git )]/ = ^(tt ^ 1 + n~J- +1 ) > 0; and b it = g ijt -igi t TT ~ t 1 > 0. Ache that A i3 

contains only one entry un when J = 2. 

As the Fisher information matrix, F is always positive semi-definite 
and thus |F| > 0 (Fedorov, 1972). As a special case, A* is the Fisher 
information at the experimental setting x* (also known as a design point or 
support point ) and thus positive semi-definite too. 

2.2 Determinant of Fisher information matrix 

Among different criteria for optimal designs, D-criterion looks for the 
allocation maximizing |F|, the determinant of F. Here, a D-optimal design 
with m predetermined design points xi,..., x. m could either be an integer¬ 
valued allocation (ni, ri 2 , ..., n m ) maximizing |F| with hxed n = n * > 
0, known as an exact design ; or a real-valued allocation (pi,p 2 , • • • ,Pm) 
maximizing |n - 1 F| with yy = n*/n > 0 and = 1, known as an 

approximate design. 

Theorem 2.2. The determinant of the Fisher information matrix 



F 





CtlH-h Ot-m — d-\-J— 1 


n. 


• • • n 


(Am 

m 
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which is an order-(d + J — 1) homogeneous polynomial of (ni,..., n m ) and 

= Y, l A ^l ( 2 -3) 

The proof of Theorem 2.2 is relegated into the Supplementary Materials 
(Section S.3). Given a map r : {1, 2,..., d + J — 1} —* {1,..., m}, A r in 
(2.3) is a (d + J — 1) x (d + J — 1) matrix whose kth row is the same as 
the kth row of A r ( fc ) , k — 1,..., d + J — 1. We denote r G (ay,..., a m ) if 
ctj = : r(j) = i} for each i — 1,..., m. 

In order to obtain analytic properties of |F|, we need a few lemmas 
including the following one, which covers Lemma 1 in Perevozskaya et al. 
(2003) as a special case: 

Lemma 2.1. Rank (A,;) = Rank(Aj 3 ) = J—1. Furthermore, is positive 
definite and 

j -i j 

i A d=n hi n ^ 1 > 0 

S=1 t=l 

where g is = (fi' _1 )'(6» s - xf/3) > 0 for s = 1,..., J - 1. 

Example 2.3. Suppose d = 2, J = 3 with link function g. According 
to Theorem 2.2, |F| in this case is an order-4 homogeneous polynomial of 
(ni,..., n m ). Based on Lemma S.4 and Lemma S.5 in the Supplementary 
Materials (Section S.2), we can remove all the terms in the form of nf, nfrij, 
or nfn 2 j from |F|. Therefore, 

m 

| F | ''f ^ ^ Cijk * U/ flj lt-h- T ''f ^ Cijki * / 1pi pl /,■ /11 

i= 1 j<k,j^i,kj^i i<j<k<l 

for some coefficients cy,-*, and Cijki ■ □ 

Based on Lemma S.4 and Lemma S.5, in order to keep 0, 

the largest possible ccj is J — 1 and the fewest possible number of positive 
ctj’s is d + 1. As a direct conclusion of Lemma S.5, the following theorem 
states that a minimally supported design has at least d + 1 support point, 
which could be much less than the number of parameters d+ J — 1: 

Theorem 2.3. In order to keep the Fisher information matrix positive 
definite, we need at least d +1 experimental settings. In other words, |F| > 0 
only if m> d + 1. 
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In order to find out whether d + 1 experimental settings or support 
points are enough to keep the Fisher information matrix positive definite, 
we study the leading term of |F| with maxi<j< m aj — J — 1. For example, 
a io = J—1 for some 1 < i 0 < m. Due to Lemma S.5 and Y^i=i a i = d+J—1, 
in order to keep c ai) ... )Q , m ^ 0, there must exist 1 < ii < 12 < • • • < id < wn 
which are different from i 0 , such that, — ■ ■ ■ — cti d = 1. The following 
lemma provides an explicit formula for such a coefficient c aij .„ )am : 

Lemma 2.2. Suppose a io = J — 1 and a ix — ■ ■ ■ — a id — 1. Then 


d 

Cai,...,am 11 &is ’ |-^-io31 ' 1-^-1 [ffi; ^1 j • • • j D] I 
s=l 

where Xi = (1 X) is an mx (d+ 1) matrix with 1 = (1,..., 1) T , X = (xi, 
..., x m ) T , and Xi [io, ii, ..., id] is the sub-matrix consisting of the ioth, Lth, 
..., irfth rows of X x . 

The proof for Lemma 2.2 is relegated to the Supplementary Materi¬ 
als (Section S.3). For the purpose of finding D-optimal allocations, we 
write |F| = /(n x , ..., n m ) for an order-(d + J — 1) homogeneous poly¬ 
nomial function /. The D-optimal exact design problem is to find an 
integer-valued allocation (n x ,..., n m ) maximizing /(n x ,..., n m ) subject to 
n l G {0,1,..., n}, i — 1,..., m and n x + • • • + n m = n with given positive 
integer n. Denote pi = ni/n, i = 1,... pm. According to Theorem 2.1, 


f(ni, ■ ■ -,n m ) 


m 

i= 1 


m 

n Pi A, 

i =1 


= n d+J - 1 f(p 1 ,...,p m ) (2.4) 


Due to (2.4), Theorems 2.2 and 2.3 can be directly applied to approximate 
designs too, while the D-optimal approximate design problem is to find a 
real-valued allocation (p x ,... ,p m ) maximizing f(pi,P 2 , ■ ■ ■ ,Pm ) subject to 
0 < Pi < 1, i — 1,..., m and p x + ■ —I- p m — 1. 

According to Lemma 2.1, |A* 0 3 1 > 0. Thus c aii ... iam in Lemma 2.2 is 
positive as long as X x [i 0 , ■ ■ ■ ,id\ is of full rank. Theorem 2.3 implies that 
a minimally supported design contains at least d + 1 support points, while 
the following theorem states a necessary and sufficient condition for the 
minimum number of support points to be exactly d + 1: 

Theorem 2.4. The Fisher information matrix is positive definite for some 
approximate design if and only if the extended design matrix X x = (1 X) 
is of full rank. In other words, /(p) > 0 for some p = (p±,... ,p m ) T if and 
only if Rank(Xi) — d+ 1. 
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The result stated in Theorem 2.4 appears to be mysterious at the first 
sight, since it claims that the minimal number of experimental settings re¬ 
quired could be strictly less than the number of parameters. We use the 
motivating example presented in Table 1 to interpret the reason. In the 
odor removal study, the main-effects cumulative link model (2.1) involves 
four independent parameters - two /Ts for the covariates (d = 2) and two d’s 
for the intercepts (J — 1 = 2). According to Theorem 2.4, a minimally sup¬ 
ported design may involve only three experimental settings. This apparent 
anomaly is due to the fact that, for multinomial responses with J = 3 cate¬ 
gories, we get two degrees of freedom (i.e., Yn,Y i2 ) from each experimental 
setting. Due to the same reason, the optimal allocation of experimental 
units is often not uniform (see Section 4), contrary to the case of binary 
responses (Yang, Mandal, and Majumdar, 2016; Yang and Mandal, 2015). 

3. D-optimal Approximate Design 

A (locally) D-optimal approximate design is a real-valued allocation p = 
{p i,... ,p m ) T maximizing /(p) = f(pi ,... ,p m ) with pre-specified values of pa¬ 
rameters. The solution always exists since / is continuous and the set of feasible 
allocations S := {(pi, ... ,p m ) T € M' m | p t > 0 ,i = 1, YiLiPi = 1} is 

convex and compact. Theorem 2.4 ascertains that a nontrivial D-optimal ap¬ 
proximate design problem needs the following assumption: 

Assumption 3.3. There are at least d+1 experimental settings and the extended 
design matrix Xi is of full rank. In other words, m > d+1 and Rank(Xi) = d+1. 

Assumption 3.3 is assumed for the rest of the paper. Then the set of valid 
allocations S + := {p = (pi,... ,p m ) T £ S | f(p) > 0} is nonempty. Since 
F = Yl=i n iA-i = n YaL\ PiA-i is linear in p and </(•) = log j - | is concave on 
positive semi-definite matrices, then /(p) = n 1 ~ d ~ J |F| is log-concave (Silvey, 
1980) and thus S+ is also convex. The following theorem characterizes all valid 
allocations in terms of the extended design matrix Xi: 

Theorem 3.5. A feasible allocation p = (pi,... ,p m ) T satisfies /(p) > 0 if and 
only if Rank(Xi [{/ | jy > 0}]) = d + 1, where Xi[{* | p % > 0}] is the sub-matrix 
consisting of the {i \ p, > 0 }th rows of Xi . In other words, S + = { p = 
(pi,P2, ■ ■ ■ j Pm) T € S | Rank(Xi[{* | p, > 0}]) = d + 1}. 

As a direct conclusion of Theorem 3.5, S+ contains all p whose coordi¬ 
nates are all strictly positive. A special case is the uniform allocation p n = 
(1/m,..., l/m) T . 

A necessary and sufficient condition for an approximate design to be D- 
optimal is of the general-equivalence-theorem type (Kiefer, 1974; Pukelsheim, 
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1993; Atkinson et al., 2007; Stufken and Yang, 2012; Fedorov and Leonov, 2014; 
Yang, Mandal, and Majumdar, 2016), which is convenient when searching for 
numerical solutions. Following Yang, Mandal, and Majumdar (2016), for a given 
P = (.Ph ■ ■ ■, Pm) T £ S. |_ and i £ {1,..., m}, we define 


fi(z ) = / 


Z -Pi,---,T — —Pi-i, z,t — —Pi+ 1 > 


1 ~Pi 


1 ~Pi 


1 ~Pi 


1 - z 
1 - Pi 


Pm 


(3.1) 


with 0 < z < 1. Note that fi(z) is well defined as long as pi < 1. Following the 
proof of Theorem 2.4, we obtain the next theorem on the coefficients of fi(z): 

Theorem 3.6. Suppose p = (pi,... ,p m ) T £ S + and i £ {1,..., m}. For 0 < 
z < 1, 

J -1 

m = (i - i - (3-2) 

i=o 

w/iere a 0 = /*(0), (aj_i,... ,ai) T = Bj^c, Bj_i = (s t_1 ) st is a (J-l) x (J-l) 
matrix, and c = (ci,..., cj_i) T withcj = (j+l) d+J ~ 1 j~ d fi(l/(j+l))-j J ^ 1 fi(0). 


Following the lift-one algorithm proposed in Yang, Mandal, and Majumdar 
(2016), we have parallel results and algorithm for our case. For simplicity, we 
also call our algorithm lift-one algorithm. 

Theorem 3.7. Given an allocation p = (p] ,... ,p^f) T € S+ , p is D-optimal if 
and only if for each i = 1,..., m, fi(z), 0 < z < 1 attains it maximum at z = p* . 


Our lift-one algorithm: 

1° Start with an allocation Po = (pit ■ ■ ■ ,Pm) T satisfying / (po) > 0. 


2° Set up a random order of i going through {1,2,..., m}. 


3° For each i, determine fi(z) according to Theorem 3.6 with J determinants 
fi( 0), fi( 1/2), /*(l/3), ..., /i(l/J) calculated according to (3.1). 

4° Use quasi-Newton method with gradient defined in (S.13) to find z* max¬ 
imizing fi(z ) with 0 < z < 1. If fi(z*) < fi( 0), let z* = 0. Define 
pi t} = {pi(l - z*)/( 1 - pi), Pi- i(l - z*)/( 1 - Pi), z*,p i+ i(l - z*)/(l - 
Pi), ■ ■ ■,Pm (1 - z*)/(l ~Pi)) T ■ Note that /(pl l) ) = fi{z*)- 

5° Replace po with p*\ and / (po) with /(p*^). 

6° Repeat 2° ~ 5° until convergence, that is, /(po) = /(pi*'*) for each i. 


Theorem 3.8. When the lift-one algorithm converges, the resulting p maximizes 
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Example 3.4. Odor removal study In our motivating example (Section 1), 
the response was ordinal in nature, serious odor, medium odor, and no odor. 
We consider the logit link and fit the cumulative link model (2.1) to the data 
presented in Table 1. The estimated values of the model parameters are (/%, 
/?2 , Oi, # 2 ) T = (—2.44,1.09, —2.67, —0.21) 2 . Suppose a follow-up experiment 
is planned and the estimated parameter values are regarded as the true value. 
The D-optimal approximate allocation found by the lift-one algorithm is p c = 
(0.4449,0.2871, 0, 0.2680) r . The efficiency of the uniform one p„ = (1/4,1/4, 
1/4, 1/4) T is (/(p«)//(p 0 )) 1 ^ 4 = 79.7% which is far from satisfactory. □ 

Example 3.5. Wine bitterness study Christensen (2015, Table 1) aggre¬ 
gated the wine data from Randall (1989). It contains the output of a facto¬ 
rial experiment with two treatment factors each at two levels (Temperature x \: 
cold (—) or warm (+); Contact X 2 ■ no (—) or yes (+)) affecting wine bitter¬ 
ness. The response was ordinal with five levels (from “1” being least bitter 
to “5” being most bitter). The original design employed a uniform allocation 
P« = (1/4,1/4,1/4,1/4) The estimated parameter values under the logit link 
are (/3i, /? 2 , 0i, @2, $ 3 , 04 ) T = (1.25,0.76, —3.36, —0.76,1.45, 2.99) T . If a follow-up 
experiment is planned regarding the estimated values of the parameters as the 
true values, then the D-optimal approximate allocation found by the lift-one al¬ 
gorithm is p 0 = (0.2694, 0.2643, 0.2333, 0.2330) t . The efficiency of the original 
design p„ is 99.9%. Nevertheless, the corresponding efficiency may drop to 80% if 
|/5i| and |/? 2 1 are both larger than 3 (see Figure 1(a)). In that case, the D-optimal 
allocations are minimally supported, that is, some pi is 0 (see Figure 1(b)), which 
will be discussed further in Section 5. □ 

With all examples that we study, the lift-one algorithm often converges 
within a few iterations. Yang, Mandal, and Majumdar (2016) also provided 
a modified lift-one algorithm, which is guaranteed to converge. The same proce¬ 
dure could be easily applied to the lift-one algorithm here if it does not converge 
in a pre-specified number of iterations. 

4. D-optimal Exact Design 

A (locally) D-optimal exact design is an integer-valued allocation n = (m, 

..., n m ) T maximizing |F| with given parameter values and a pre-specified total 
number n of experimental units or runs. In the design literature, different dis¬ 
cretization methods have been proposed to round an approximate design into an 
exact design for a given n, including the quota method (Kiefer, 1971; Pukelsheim, 
1993) and the efficient rounding procedure (Pukelsheim, 1993; Pukelsheim and 
Rieder, 1992), which usually work well for large enough n but with no guarantee 
for small sample size (Irnhof, Lopez-Fidalgo, and Wong, 2001). 
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(a) 


(b) 




P, 


P, 


Figure 1: Wine bitterness study with assumed true parameter values (/3 i,/ 32, —3.36, 
—0.76, 1.45,2.99) t : (a) contour plot of efficiency of the original design; (b) regions for a 
D-optimal design to be minimally supported, that is, pi = 0 for i = 1, 2, 3, or 4. 

In this section, we aim to provide a direct search for D-optimal exact designs. 
As a conclusion of Theorem 3.5, the following corollary implies that we need at 
least d + 1 experimental units to make |F| > 0 possible: 

Corollary 4.1. |F| > 0 if and only if RankfXi [{< | m > 0}]) = d+ 1. 

In order to avoid trivial cases, we assume n > d + 1 throughout this sec¬ 
tion. To maximize /(n) = /(m,..., n m ) = |F|, we adopt the idea of exchange 
algorithm which was first suggested by Fedorov (1972). The exchange algorithm 
used here is to adjust n,; and rij simultaneously for randomly chosen (i,j) while 
keeping n* + rij = c as a constant. 

We start with an n = (m,... ,n m ) T satisfying /(n) > 0. Following Yang, 
Mandal, and Majumdar (2016), for 1 < i < j < m, we dehne 

fij(z ) = / (ni,... ,ni-i,z,n i+ i ,.. .,rij-i,c- z,n j+1 ,... ,n m ) (4.1) 

where c = n* + n 3 , z = 0,1,..., c. Note that = /(n). From Theorem 2.2, 

Lemmas S.4 and S.5, we obtain the following formula on calculating fij(z): 

Theorem 4.9. Suppose n = (m,... ,n m ) T satisfies /(n) > 0 and ni + nj > J 
for given 1 < i < j < m. For z = 0,1,..., n* + rij , 

J 

fij(z) = ^2°sZ s 

s=0 


(4.2) 
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where cq = (0), and c±,... ,cj can be obtained using (ci,..., cj) 7 = B^di, 

..., dj) r with Bj = ( s t ~ 1 ) s t as a J x J matrix and d s = (fij(s) — /ij(0))/s. 

Note that the J x J matrix Bj in Theorem 4.9 shares the same form of 
in Theorem 3.6. According to Theorem 4.9, in other to maximize fij(z) 
with z = 0, 1 ,..., rii + nj , one can obtain the exact polynomial form of fij(z) 
by calculating fij(0), fij( 1),... , fij(J)- There is no practical need to find out the 
exact form of fijiz ) if n % + nj < J since one may simply calculate fij(z) for each 
z. Following Yang, Mandal, and Majumdar (2016), an exchange algorithm (see 
the Supplementary Materials, Section S.5) based on Theorem 4.9 could be used 
to search for a D-optimal exact allocation. 

Example 3.4 Odor removal study ( continued ) Suppose we want to con¬ 
duct a follow-up experiment with n experimental units. Using the exchange 
algorithm mentioned above, we obtain the D-optimal exact designs across dif¬ 
ferent n’s (Table 2). As expected, the D-optimal exact allocations (n-i,... ,n- 4 ) T 


Table 2: D-optimal Exact Designs and Approximate Design for Odor Removal Study 


n 

n\ 

n 2 

n 3 

714 

u _ 4 |F| 

# iterations 

Time(sec.) 

3 

1 

1 

0 

1 

0.0002911 

1 

< 0.01 

10 

4 

3 

0 

3 

0.0003133 

3 

0.02 

40 

18 

11 

0 

11 

0.0003177 

3 

0.02 

100 

44 

29 

0 

27 

0.0003180 

4 

0.05 

1000 

445 

287 

0 

268 

0.0003181 

5 

0.39 

Po 

0.4449 

0.2871 

0 

0.2680 

0.0003181 

5 

0.03 


is consistent with the D-optimal approximate allocation p G = (p\.... ,P 4 ) T (last 
row of Table 2) for large n. The time costs in seconds (last column of Ta¬ 
ble 2) are recorded on a PC with 2GHz CPU and 8GB memory. Suppose we 
rerun an experiment with n = 40. With respect to the D-optimal exact design 
n G = (18,11, 0,11) T , the efficiency of the uniform design n u = (10,10,10,10) T is 
(/(n u )//(n 0 )) 1 / 4 = 79.7% only. □ 

Example 4.6. Polysilicon deposition study Wu (2008) considered an ex¬ 
periment for studying the polysilicon deposition process with six 3-level factors, 
described in details by Phadke (1989). Due to the inconvenience of counting the 
number of surface defects, which is one of the major evaluating characteristics, 
they treated it as a 5-category ordinal variable: 1 for 0 ~ 3 defects, 2 for 4 ~ 30, 
3 for 31 ~ 300, 4 for 301 ~ 1000, and 5 for 1001 and more. The original design, 
denoted by n„ , includes 18 experimental settings based on the Lig orthogonal 
array. In order to apply a cumulative link model, we represent each 3-level factor, 
say A with levels 1,2,3, by its linear component Ai taking values —1,0,1 and 
quadratic component A 2 taking values 1,—2,1 (Wu and Hamada, 2009). Then 
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the fitted model with complementary log-log link chosen by both AIC and BIC 
criteria (see, for example, Agresti (2013)) involves 4 cut-points (tfi, a?, 03 ,dq) = 
(-1.59,-0.58,0.41,1.22) and 12 other coefficients (/?n, ^ 12 , /? 2 i> P 22 , ■ ■ ■, Pw) = 
(1.45,-0.22,1.35,0.02,-0.12,-0.34,0.19,0.00,0.22, 0.08, 0.05, 0.17). When a 
follow-up experiment is conducted and the true parameter values are assumed to 
be the estimated ones, we use the exchange algorithm and find out a D-optimal 
18-run design, denoted by n 0 (see the Supplementary Materials, Section S. 6 , for 
a list of the 18 experimental settings). Compared with n G , the efficiency of the 
original design n u is (/(n u )//(n D )) 1 / 16 = 73.1% only. In order to check the effi¬ 
ciency of a rounded design, we use the lift-one algorithm and find out that the 
D-optimal approximate design contains 100 positive pi s out of the 729 distinct 
experimental settings. In this case, both the quota method and the efficient 
rounding procedure end with the same rounded design n r (see Section S. 6 ). Not 
unexpectedly, its efficiency is (/(n r )//(n 0 )) 1//16 = 86 . 1 % only. □ 

5. Minimally Supported Design 

A minimally supported design is a design with the minimal number of sup¬ 
port points (i.e., x* with pi > 0 or n* > 0) while keeping the Fisher information 
matrix positive definite, that is, |F| >0. It is of practical significance since the 
experiment runs with the minimal number of different settings. For example, the 
18 experimental settings in the polysilicon deposition study (Example 4.6) had 
to be run in a sequential way and only two settings were arranged in each day 
(Phadke, 1989). It is usually relatively easier to repeat the experiments with the 
same setting. Less experimental settings often indicate less time and less cost. 
Another practical application of a minimally supported design is to serve as the 
initial design for numerical searching algorithms. Since the objective function in 
this case takes a much simpler form, an optimal allocation restricted to those 
support points can be obtained more easily or even analytically. 

According to Theorem 2.3, a minimally supported design contains at least d+ 
1 support points. On the other hand, according to Theorem 3.5 and Corollary 4.1, 
a minimally supported design could contain exactly d + 1 support points if the 
extended design matrix Xi = (1 X) is of full rank, that is, Rank(Xi) = d + 1. 

Example 5.7. Let J = 2 and then the multinomial response is actually binomial. 
In this case, there are d+ 1 parameters, 0i,/?i,... ,f3d . For a general link function 
g satisfying Assumptions 2.1 and 2 . 2 , gi 0 = gi 2 = 0 , gn = (g~ l )'(9 1 — xf(3) > 0 , 
e* = un = cn = g‘f 1 /['n i i(l - 7r*i)], i = l,...,m. Then A i3 in Theorem 2.1 
contains only one entry, m 1 , and thus | A* 3 1 = u t 1 or simply e t (Lemma 2.1 still 
holds). Assume further that the mx d design matrix X satisfies Assumption 3.3. 
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According to Theorem 2.2, Lemmas S.4, S.5, and 2.2, given p = (pi, ■ ■ ■ ,p m ) T , 
/(P) = n' (d+1) |F| = |Xi[* 0 ,n,.. ■ iidlfpioeioPhCii ''' Pi d e i d 

l<io<ii<—<i d <m 

(5.1) 

Note that (5.1) is essentially the same as Lemma 3.1 in Yang and Mandal (2015). 
Then a minimally supported design may contain d + 1 support points and a D- 
optimal one keeps equal weight l/(d+l) on all support points (Yang and Mandal, 
2015, Theorem 3.2). □ 

For univariate responses (including binomial ones) under a generalized linear 
model, a minimally supported design must keep equal weights on all its support 
points in order to keep D-optimality (Yang and Mandal, 2015). However, for 
multinomial responses with J > 3, it is usually not the case. In this section, we 
use one-predictor (d = 1) and two-predictor (d = 2) cases for illustration. 

In order to check if a minimally supported design is D-optimal, we need 
a Karush-Kuhn-Tucker-type condition. Since /(p) is log-concave, the Karush- 
Kuhn-Tucker conditions (Karush, 1939; Kuhn and Tucker, 1951) are also suffi¬ 
cient. We have the following theorem as a direction conclusion. 

Theorem 5.10. An allocation p = (pf,... ,p^) T satisfying /(p) > 0 is D- 
optimal if and only if there exists a A G R such that df(p)/dpi = A if p* > 0 or 
< A if p* = 0, i = 1,..., m. 

5.1 Minimally supported designs with one predictor 

We first consider the cases with one predictor (d = 1). We assume the 
response has at least three categories (J > 3) to skip the well-known cases 
(Example 5.7). The corresponding parameters here are /3i and 9 1 ,..., 0j_i . 

We start with designs supported on two points (rn = 2, minimally sup¬ 
ported). As a direct conclusion from Theorem 2.2, Lemmas S.4 and S.5, we 
have: 

Theorem 5.11. Suppose d = 1, J > 3, m = 2. The objective function 

J-i 

f{Pi,P 2 ) = n' 2 |F| = Y c sp{~ s pI (5.2) 

S= 1 

where (ci,..., cj_ i) T = • • •, dj-\) T with B./_i = ( s t ~ 1 ) s t as a (J — 1) x 

(J — 1) matrix and d s = /(l/(s + 1), s/(s + 1)) • (s + 1) J /s. 

Actually, according to Lemma 2.2, c\ = e 2 n s=i 9i s ■ nil 7r it'( x 1 ®2) 2 , 
Cj-i = e\ nii^-nii n 2t 1 (x 1 — X 2 ) 2 , where X\,X 2 are the predictor levels. 
Theorem 5.11 provides a practically useful way to find out the exact form of 
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f(PhP 2 ) after calculating |F| for J — 1 different allocations. Then the D-optimal 
problem is to maximize an order-J polynomial f(z, 1 — z) for z € [0,1]. As a 
special case, the D-optimal allocation of J = 3 can be solved explicitly as follows: 

Corollary 5.2. Suppose d = 1, J = 3, m = 2. The objective function 


f(pi,P 2 ) = PlP 2 (ciPl + C 2 P 2 ) (5.3) 

where c L = e2gi l gi 2 {^u^i2^i3)~ 1 (xi - x 2 ) 2 > 0 , c 2 = eis^s^^iTr^Tr ^) -1 
(xi — x 2 ) 2 > 0, and xi,x 2 are f/ie two levels of the predictor. The D-optimal 
design p = (p{,p 2 ) 

Cl - C 2 + y/cf - ClC 2 + C 2 Cl 

Pi — “ / o 9 > P2 — - / 9 o 

2Ci - c 2 + A/cf - ClC 2 + C 2 2cj. - c 2 + y/cf ~ ClC 2 + c 2 

Furthermore, p\ = p 2 = 1/2 */ and only if c\ = c 2 . 

Under the setup of Corollary 5.2, it can be verified that p\ = p 2 = 1/2 if 
/3i = 0. In general p* / p 2 , and p\ > p 2 if and only if ci > c 2 . The following 
result provides conditions for D-optimality of such a minimally supported design. 
Its proof is relegated to the Supplementary Materials (Section S.3). 

Corollary 5.3. Suppose d = 1, J = 3, m > 3. Let xi,... ,x m be the m distinct 
levels of the predictor. A minimally supported design p = (p*,p 2 ,0,..., 0) T is 
D-optimal if and only if 

(1) p\,p* 2 are defined as in (5.f), 

(2) saiPi) 2 + (si 5 - 2ci)p|p^ + (s i4 - c 2 ){p* 2 ) 2 < 0, i = 3,... ,m, 

where ci,c 2 are same as in Corollary 5.2, 5^3 = ej< 7 2 1 p 2 2 ( 7 rii 7 ri 2 7 ri 3 )' 1 (xi — 
Xj ) 2 > 0 , s i4 = ei< 7 21 < 7 2 2 ( 7 r 2 i 7 r 22 7 r 23 )~ 1 (x 2 - x *) 2 > 0 , s i5 = ei(u 2 2 Un + U 2 iu i2 - 
2622 ^ 2)^1 - a: 2 )(xi - Xj) + e 2 (ui 2 Un + ttntt j2 - 26i 2 6j 2 )(x 2 - xi)(x 2 - xf) + 
ei(u 12 U 21 + U U U22 - 2b V2 b22)(Xi - Xi)(Xi - x 2 ). 

Example 5.8. Suppose d = 1, J = 3, and m = 3 with factor levels {—1,0,1}. 
Under the logit link g, the parameters (3,9\, 0 2 satisfy gC/\j) = 0j+P, g(l 2 j) = Qj, 
g(C'ij) = @j — P, j = 1,2. We investigate when a D-optimal design is minimally 
supported, that is, Pi = 0 for some i. According to Theorem 5.11, a D-optimal 
deign satisfies pi = p 3 = 1/2 if P = 0. Figure 2 shows cases with more general 
parameter values. In Figure 2(a), four regions in (9\, 6 ) 2 )-plane are occupied by 
minimally supported designs (9\ < 0 2 is required). For example, regions labeled 
with p 2 = 0 indicates a minimally supported design satisfying p 2 = 0 is D-optimal 
given such a triple {9\,02,P = —2). From Figure 2 (b), a design supported on 
{—1,1} (that is, P 2 = 0) is D-optimal if P is not far away from 0. □ 
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(a) 0=—2 


(b) e 2 =5 




P 


Figure 2: Regions for a two-point design to be D-optimal with d = 1, J = 3, x G 
{ — 1,0,1}, and logit link (note that Q\ < 82 is required) 


Example 5.9. Toxicity study Agresti (2013, Table 8.7) reported data from a 
developmental toxicity study with one factor (concentration of diEGdiME at five 
levels: 0, 62.5, 125, 250, 500 mg/kg per day) and a 3-category ordinal response 
(status of mouse fetus: nonlive, malformation, or normal). In this case, d = 1, 
J = 3, and m = 5. For illustration purpose, we fit a cumulative link model 
with cauchit link chosen by both AIC and BIC criteria. The estimated pa¬ 
rameter values are (/?i,0i,02) 7 = (—0.0176,-8.80,-5.34). Suppose a follow-up 
experiment is planned with (/%,#!, # 2 ) r regarded as the true parameter values, 
then the D-optimal approximate allocation found by the lift-one algorithm is 
p 0 = (0,0, 0,0 .4285,0.5715) r , which is minimally supported. Alternatively, for 
each pair of indices 1 < i < j < 5, we obtain the best design ( p*,pj ) 

supported only on Xi,Xj according to Corollary 5.2, then check whether (p* , p *) 
is D-optimal using Corollary 5.3. We confirm that p D is the only minimally sup¬ 
ported design that is also D-optimal. With respect to p D , the efficiency of the 
original design (roughly a uniform one) is only 52.6%. □ 

5.2 Minimally supported designs with two predictors 

In this section, we consider experiments with two predictors (d = 2) and a 
three-category response (J = 3). The parameters are Pi, fh, 8 \, 82 . For cases 
with J > 4, similar conclusions could be obtained but with messier notations. 

According to Theorem 2.3, a minimally supported design in this case needs 
three support points, for example, (xn,Xi 2 ), i = 1,2,3. Under Assumption 3.3, 
the 3x3 matrix Xi = (1 X) is of full rank. Following Theorem 2.2, Lemmas S.4, 
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S.5, and 2.2, the objective function with ( d,J,m ) = (2,3,3) is 

f(PuP2,P3) = |Xi| 2 eie 2 e 3 ■ PiP2P3(wiPi + W 2 P 2 + w 3 p 3 ) (5.5) 

where w t = e ~ 1 </ 2 ,1+27G.3)~ 1 > 0. Since f(p 3 , p 2 , p 3 ) = 0 if P 1 P 2 P 3 = 0, we 
only need to consider p = ( Pi,P 2 ,P 3 ) T satisfying 0 < Pi,P 2 ,P 3 < 1- 
According to Theorem 5.10, p maximizes f(p \, p 2 ,p 3 ) only if 

di = di = di 

dpi dp 2 d P3 [ ' 1 

Following Tong, Volkmer, and Yang (2014), we obtain its analytic solution: 

Theorem 5.12. Without loss of generality, we assume w\ > W 2 > w 3 > 0 . The 
allocation p = (Pi,P 2 ,P 3 ) T maximizing ,f(pi-/P 2 ,P:i) in (5.5) exists and is unique. 
It satisfies 0 < p% < p\ < p\ < 1 and can he obtained analytically as follows 

(i) If wi > iu 2 = w 3 , then p\ = Ai/(4wi + A 3 ), p* 2 = p £ = 2wi/{Awi + Ai), 
where A\ = 2w 3 — 3 w 2 + \/ 4 w\ — A.W 1 W 2 + 9w 2 • A special case is p\ 
= P 2 =P *3 = V 3 */ wi=w 2 = w 3 . 

(ii) Ifwi = W 2 > w 3 , then p\= p* 2 = A 2 /[2(A 2 + 2wi)\, p 3 = 2 ioi/(A 2 + 2wi), 
where A 2 = 3w\ — 2 w 3 + y / 9rcf — Aw\w 3 + 4tc 3 . 

(in) If wi > w 2 > w 3 , then p\ = yi/(yi + y 2 + 1), p\ = 2/2/(271 + 2/2 + 1), 
P* 3 = 1/(2/1 + 2/2 + 1), where 

_ b 2 2 1 / 3 (3&i - b|) A 1/3 _ (wi - w 3 )j/i 

^ 3 3A 1 / 3 3 x 2 1 / 3 ’ (w 2 — w 3 ) + (u)] — w 2 )yi 

with A = —27bo + 9b\b 2 — 26 3 + 3 3 / 2 (27&q + 46 3 — I 860&1&2 — fr 2 6 2 + 46o& 2 ) 1//2 , 
bi = Ci/c 3 , i = 0,1,2, and c 0 = w 3 (w 2 - w 3 ) > 0, c\ = 3wiw 2 - w 3 w 3 - 
Aw 2 w 3 + 2w‘l > 0, c 2 = 2w\ — Aw\W 2 — wiw 3 + 3w 2 w 3 , c 3 = w\{w 2 — wi) < 0. 

The proof of Theorem 5.12 is relegated to the Supplementary Materials (Sec¬ 
tion S.3). Again, a minimally supported design is usually not uniformly allocated. 

Corollary 5.4. Suppose d = 2, J = 3, and m = 3. Then p = (1/3,1/3, l/3) T is 
D-optimal if and only if w 1 = w 2 = w 3 , where w\,W 2 ,w 3 are defined as in (5.5). 

Example 5.10. Consider a typical 2 2 factorial design problem with a three- 
category response and four design points (1,1), (1,-1), (—1,1), (—1,-1), denoted 
by (, xn,Xi 2 ),i = 1,2,3,4. That is, d = 2, J = 3, and m = 4. Define Wi = 
er^iSafaiWa) -1 , * = 1,2, 3,4. There are five special cases: (i) if fii = = 

0, then wi = W 2 = w 3 = W 4 ; /ii/ if = 0,/? 2 / 0, then rci = rc 3 , rc 2 = 104, but 
/ rc 2 ; /iiij if fii / 0, /3 2 = 0, then tci = w 2 , tc 3 = 104, but rci / u; 3 ; (iv) if 
/3i = /+2 / 0, then u; 2 = rc 3 , but w\,w 2 ,W 4 are distinct; (v) if /3i = —/3 2 / 0, 
then wi = 1U4, but w\,W 2 ,w 3 are distinct. □ 
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Theorem 5.12 provides analytic forms of minimally supported designs with 
d = 2 and J = 3. As a direct conclusion of Theorem 5.10, we have the necessary 
and sufficient conditions for D-optimality of a minimally supported design: 

Corollary 5.5. Suppose d = 2, J = 3, and m > 4. Let (xju x^), i = 1,..., m be 
rn distinct level combinations of the two predictors. Then Xi = (1 X) is anm,x3 
matrix. A minimally supported design p = (Pi,p 2 ,p 3 , 0, ..., 0) T is D-optimal if 
and only if (1) p*, p 2 , p 3 are obtained according to Theorem 5.12; and (2) 

|Xi [1, 2, i\ \ 2 e\e2eip\p2(wip\ + w 2 p* 2 ) + \^i[l,3,i]\ 2 eiezeip\p%{wip\ + w 3 p^) 
+ |Xi [2, 3, i]\ 2 e 2 e 3 eip* 2 p* 3 {w 2 p* 2 + w 3 p 3 ) + A/hICJH 
< |Xi[l, 2 , 3 ]| 2 eie 2 e 3 p 2 P 3 ( 2 u, iPi+ W 2 P 2 + ^ 3 ^ 3)1 for i = 4,... ,m, 

where ej = u 3 \ + u j2 - 2 b j2 , w 3 = g 2 jX g 2 j 2 {jr 32 e 33 )~ x , j = 1 Di = 

j,k, s,t}£Ei e j e k{u s iu t2 + u s2 uti - 2b s 2 b t2 ) • |Xi[j, k,s\ \ ■ |Xi[j, k,t\\ with the sum 
over Ei = {(1,2,3,*), (1,3,2,*), (1,*,2,3), (2,3,1,*), (2,*, 1,3), (3,*, 1,2)}. 

Example 5.11. Consider experiments with two predictors, a three-category re¬ 
sponse, and four design points, modeled with logit link. That is, d = 2, J = 3, 
rn = 4. Suppose the four design points are (1, 1), (1, —1), (—1, 1), (—1,-1). 
According to Theorem 5.12 and Corollary 5.5, we can analytically calculate the 
best minimally supported design (that is, supported on three points) and deter¬ 
mine whether it is D-optimal or not. Figure 3 provides the boundary lines of 
regions of parameters (fi\,fi 2 .0 \, 0 2 ) for which the best three-point design is D- 
optimal. In particular, Figure 3(a) shows the region of (/3i,f3 2 ) for given 0\. 0 2 ■ 
It clearly indicates that the best three-point design tends to be D-optimal when 
the absolute values of /3 i,/ 3 2 are large. The region tends to be larger as the 
absolute values of 6 i, 6 2 increase. On the other hand, Figure 3(b) displays the 
region of (Q\ ,d 2 ) for given (A , fi 2 . The symmetry of the boundary lines about 
0\ + 9 2 = 0 is due to the logit link which is symmetric about 0. An interesting 
conclusion based on Corollary 5.5 is that in this case a three-point design can 
never be D-optimal if = 0 or f3 2 = 0. □ 

Remark 5.1. In this section we show that the extra degrees of freedom plays 
an important role against the uniformity of D-optimal allocation in a minimally 
supported design. For multinomial-type responses with J categories, the total 
degrees of freedom from m distinct experimental settings is m(J — 1), while a 
cumulative link model contains d + J — 1 parameters. For a minimally supported 
design, rn = d+1 and m(J— 1) = d+J — 1 if and only if J = 2 (see Example 5.7). 
In that case, the objective function /(p) oc Pi 0 p %1 ■ ■ - Pi d and thus the D-optimal 
allocation is pi 0 = Pn = • ■ ■ = Pi d = l/(d+l). However, if J > 3, then the degrees 
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(a) 


(b) 




Figure 3: Boundary lines for a three-point design to be D-optimal with logit link: Region 
of (/3i, fa) for given (0i, 62) is outside the boundary lines in Panel (a); Region of (0i, 62) 
(with 9\ < 62 ) for given (/3i,/32) is between the boundary lines and 6 \ = 62 in Panel (b) 

of freedom is strictly bigger than the number of parameters and there are “extra” 
degrees of freedom. In this case, distinct experimental settings may play different 
roles in estimating the parameters values. For example, if d = 1, J = 3, m = 2, 
the objective function /(p) = piP 2 (c\p\ + C2P2) according to Corollary 5.2; if d = 
2, J = 3, m = 3, /(p) oc p\p 2 P‘i(wiPi + W 2 P 2 + W 3 P 3 ) according to equation (5.5). 
In those cases, the D-optimality of a uniform allocation depends on ci = C 2 or 
w\ = W 2 = W 3 , which is not true in general. 

6. EW D-optimal Design 

The previous sections mainly focus on locally D-optimal designs which re¬ 
quire assumed parameter values, (/?i,..., f3d, 9i, ... ,9 j_i). For many applica¬ 
tions, the experimenter may have little or limited information about the values 
of parameters. In that case, Bayes D-optimality (Chaloner and Verdinelli, 1995) 
which maximizes E(log |F|) given a prior distribution on parameters provides a 
reasonable solution. An alternative to Bayes one is the EW D-optimality (Yang, 
Mandal, and Majumdar, 2016; Atkinson et ah, 2007) which maximizes log |E(F)| 
essentially. According to Yang, Mandal, and Majumdar (2016)’s thorough simu¬ 
lation study across different models and choices of priors, EW D-optimal designs 
are much easier to calculate and still highly efficient compared with Bayes ones. 

Based on Theorem 2.1, an EW D-optimal design which maximizes |E(F)| 
may be viewed as a locally D-optimal design with ei,Cit,uu and bn replaced by 
their expectations. After the replacement, Lemma S.2 still holds. Therefore, 
almost all the lemmas, theorems, corollaries, and algorithms in the previous 
sections can be applied directly to EW D-optimal designs as well. The only 
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exception is due to Lemma 2.1 which provides the formula of A j 3 1 in terms of 
gij and nij . In order to fit EW D-optimal designs, |A i3 | needs to be calculated 
in terms of un and bn . For example, |A,; 3 | = un if J = 2, |A,; 3 | = unUj .2 — &f 2 ^ 
J = 3, and |A,; 3 | = unu^Uis — unb? 3 — tii 3 6 f 2 if J = 4. Then the formulas of |A,; 3 | 
in Lemma 2.2, ci, c 2 in Corollary 5.2, Sj 3 , 5 ^ 4 , s *5 in Corollary 5.3, Wi in (5.5), and 
Wj in Corollary 5.5 need to be written in terms of un and bn as well. 

According to Lemma S. 2 , we only need to calculate E(un), i = 1,..., m; t = 
1,..., J — 1 and E(bn),i = 1,..., m; t = 2,..., J — 1 (if J > 3). Then E(cn) = 
E{un) — E{bit) — E(bij+i) and E(ei) = Ylt =1 After that, we can use the 

lift-one algorithm in Section 3 or the exchange algorithm in Section 4 to find EW 
D-optimal designs. We use the odor removal example to illustrate how it works. 

Example 3.4 : Odor Removal Study ( continued ) Again suppose that 
we plan to conduct a followup experiment. Instead of assuming the parameter 
values (/?i,/ 32 i^ii # 2 ) = (—2.44,1.09,-2.67,-0.21), suppose we believe that the 
true values of parameters satisfy /% G [—3, — 1], /3 2 € [0, 2], 9\ G [—4, —2], and 0 2 G 
[—1,1]. In order to perform Bayes optimality, we assume that the four parameters 
are independently and uniformly distributed within their intervals. We use R 
function constrOptim to maximize ^(p) = if(log|F|) and find the Bayes D- 
optimal allocation p^ = (0.3879,0.3264, 0.0000,0.2857) r . The procedure costs 
313 seconds computational time. In order to get the EW D-optimal design, we 
only need 5.43 seconds in total to calculate E{un ) and E(bn) and find p e = 
(0.3935,0.3259, 0, 0.2806) r using the lift-one algorithm. Even in terms of Bayes 
Optimality (Chaloner and Larntz, 1989; Song and Wong, 1998; Abebe et ah, 
2014), the relative efficiency of p e with respect to p& is exp{(<^(p e ) — 0(pb))/4} x 
100% = 99.99%, while the relative efficiency of the uniform allocation p. a = 
(0.25, 0.25,0.25,0.25) T is only 87.67%. 

In order to check robustness towards misspecified parameter values, we let 
6 = (/3i, /3 2 , 0i, 0 2 ) T run through all 0.1-grid points in [—3, —1] x [0, 2] x [—4, —2] X 
[—1,1], that is, the set of I = {—3, —2.9,..., —1} x {0,0.1,..., 2} x {—4, —3.9, 
..., —2} x {—1, —0.9,... , 1}. For each 6 G I, we use the lift-one algorithm to 
find the D-optimal allocation pe and the corresponding determinant f(pg) = 
|F(pfl)|, and then calculate the efficiency (/(p)//(pe )) 1,/4 for p = pfe,p e , and 
p u , respectively. Table 3 shows the summary statistics of the efficiencies for 
6 G I. It implies that p& and p e are comparable and both of them are much 
better than p u in terms of robustness. □ 

7. Discussion 

In this paper, we use four real experiments to illustrate how much improve¬ 
ments the experimenter could make by adopting our D-optimal designs, including 
Example 3.4: Odor removal study; Example 3.5: Wine bitterness study; Exam¬ 
ple 4.6: Polysilicon deposition study; and Example 5.9: Toxicity study. In the 
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Table 3: Summary of Efficiency in Odor Removal Study 


Design 

Min. 

1st Quartile 

Median 

Mean 

3rd Quartile 

Max. 

Bayes p b 

0.8464 

0.9813 

0.9915 

0.9839 

0.9964 

1.0000 

EW p e 

0.8465 

0.9802 

0.9917 

0.9838 

0.9967 

1.0000 

Uniform p u 

0.7423 

0.8105 

0.8622 

0.8674 

0.9249 

0.9950 


absence of theoretical results on optimal designs for ordinal data, their original 
designs are either uniform designs with equal number of replicates or regular 
fractional designs (Example 4.6) following the classical results under linear mod¬ 
els. Although these original designs are common in practice, their efficiencies 
compared with our D-optimal designs are often far from satisfactory: 79.7% in 
Example 3.4, 73.1% in Example 4.6, and 52.6% in Example 5.9. 

More interestingly, our D-optimal designs recommended for Example 3.4 
and Example 5.9 are both minimally supported, that is, these designs require 
the minimum number of experimental settings, which indicates the minimum 
cost and minimum time spent for running the experiments (Phadke, 1989). We 
have two surprising findings that are different from the cases under univariate 
generalized linear models (Yang and Mandal, 2015): (1) the minimum num¬ 
ber of experimental settings can be strictly less than the number of parame¬ 
ters: 3 experimental settings for 4 parameters in Example 3.4, and 2 experi¬ 
mental settings for 3 parameters in Example 5.9; (2) the allocation of exper¬ 
imental units on the support points of a minimally supported design is usu¬ 
ally not uniform: = (0.4449,0.2871,0.2680) in Example 3.4, and 

[PiiPb) = (0.4285,0.5715) in Example 5.9. These findings have been justified 
theoretically. The apparent anomaly is explained by the fact that a single exper¬ 
imental setting or support point with multiple-category responses generates more 
than one independent summary statistics, and thus produces enough degrees of 
freedom to estimate more parameters than the number of experimental settings. 

Even when the experimenter has little information about the parameters, we 
recommend the use of EW D-optimal designs instead of uniform designs. It is 
easy to obtain them using our algorithms and they are much more robust than 
uniform designs, which is illustrated using Example 3.4. The EW D-optimal 
designs are highly efficient surrogate to Bayesian D-optimal designs as well. 

Cumulative link models are widely used for modeling ordinal data. Nev¬ 
ertheless, there are other models used for multinomial-type responses, includ¬ 
ing baseline-category logit model for nominal response, adjacent-categories logit 
model for ordinal data, and continuation-ratio logit model for hierarchical re¬ 
sponse (see Liu and Agresti (2005), Agresti (2013) for a review). The methods 
developed in this paper could be extended for those models as well. For further 
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extension, our approaches could be used for planning experiments with more 
than one categorical responses. For example, both the paper feeder experiment 
and the PCB experiment analyzed by Joseph and Wu (2004) involved multiple 
binomial responses. The characteristics of optimal designs for these experiments 
are expected to be different from the ones with univariate responses as well. 
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Supplementary Materials 

S.l Commonly Used Link Functions for Cumulative Link Models 


Link function 
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1 — exp{—e^} 
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1 

7r( 1-P?7^ ) 


where <h _1 (-) is the cumulative distribution function of N{ 0,1), and <f{-) is the 
probability density function of N( 0,1). 

Example 2.1 ( continued ) For logit link g , g~ 1 {'i]) = e ,? /(l + e' 7 ) and ( g _1 )' = 

g~ l (l - g~ l ). Thus g tj = {g~ l )'{9j - xf/3) = (7^)(1 - 7 ij)- With J = 3, 
we have 771 + 772 + 71.^3 = 1 for i = 1 ,m. Then for i = 1,m, gn = 
77l(772 + Ki3),gi2 = (nn + VTj2)vTj3, 6j2 = (?7l + ?72)(772 + 773)>«il = 

(T;i T 772 )(t 72 T 773)7^72 — 77 "*3 7 '"j2 (^7l T 772) (772 T 773)1 Ql — 77l(77l T 
772 )(t 72 + 773), Cj2 = 77 . 3 (771 + 772) (772 + 773), Cj = (771 + 772) (771 + 773) (772 + 773) 
□ 

S.2 Additional Lemmas 

For Section 2: Since (Yu,... ,Yjj),i = 1 are independent, the log- 

likelihood function (up to a constant) of the cumulative link model is 

m J 

l{Pi, ■ ■ ■ ,Pd,9i, ■ ■ ■, 0j- 1) = EE Yij l0g(7Ty) 

i=1 3 = 1 
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where nij = 7ij “Tm-i with = g 1 (0 j -xf(3) for j = 1,..., J -1 and 7 i0 = 0, 
7 ij = 1, i = 1,..., m. For s = l,...,d, t = 1,..., J - 1, 




dl 

dd t 


m 

2=1 


— ' (5 -xf/3) 


+—•[(& -xf/3) - (3 1 ) , (0i -xf/3)] 

772 

+ • • • [-(s _ 1 )'(0J-i - xf/3)] | 


5 > -1 m 

i=l 


x f/3) 



*WA 

TTijt+l / 


Since Ty’s come from multinomial distributions, we know E(Yij ) = 77717 ,- , 
E(Y z 2 j) = m(ni - 1)7 rfj + Tii'Kij , and E(Y is Y it ) = m(ni - 1 ) 717,717 when s^t. 
Then we have the following lemma: 


Lemma S.l. Let F = (F st ) be the (d + J — 1) x (d + J — 1) Fisher information 
matrix. 


(i) For 1 < s < d, 1 < t < d, 

F = , — E ' 9 ‘ 81 


df3 s d fa 


m J / ^2 

E \ ■' \9ij ~ 9i,j-l) 

niXi s Xit ^ —---- 


i=l 


3 =1 


where gij = {g )'{0j - x/ /3) > 0 for j = 1,..., J - 1 and g i0 = g tJ = 0. 


fii) For 1 < s < d, 1 < t < J — 1, 


F s ,d+t — E 


dl dl 
d/3 s d0 t 


y^77(-37s)ff^ 


2=1 


9it 9i,t—l ^2,t+l 9it 


ft it 


fti,t-\-l 


(Hi) For l<s<J-l,l<t<d, 


Fd+s,t — E 


dl dl 
d0 s d/3 t 
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(iv) For 1 < s < J — 1, 1 < t < J — l, 

if s = t 
if\s~t\ = l 
if |s — t\ > 2 

where s V t = max{s, t}. 

Perevozskaya et al. (2003) obtained a detailed form of Fisher information 
matrix for logit link and one predictor. Our expressions here are good for fairly 
general link and d predictors. To simplify the notations, we denote 

e* = > i = (S.l) 

9it ~ 9i,t— 1 f/ijt+l fto \ 1 T 1 /Q 0 \ 

- : - : - , i = 1,..., to; t = 1,..., J — 1 (S.2) 

u it = + 7r M+i) > °> * = !,••■,»";* = if-■,«/-! (S.3) 

bu = 9i,t-i9it'Kf t 1 > 0, i = 1,... ,m; t = 2,,.., J — 1 (if J > 3) (S.4) 

Note that gij is defined in Lemma S.l (i). Then we obtain the following lemma 
which plays a key role in later on calculation of |F|. 

Lemma S.2. c it = u it -b it -bij + i, i = 1 = 1 ,..., J-l; e, = Yft=\ c n = 

(EfT, 1 (u rf - 2 b it ), i = where bn = bu = 0 for i = 1 ,..., m. 

Lemma S.3. Rank((A*i A^)) < 1 where “=” is true if and only if x, / 0. 

Based on Lemma 2.1 and Lemma S.3, we can obtain the two lemmas below 
on c ai am which significantly simplify the structure of |F| as a polynomial of 
{n \, • • •, u m ). 

Lemma S.4. If maxi> J, then |A r | = 0 for any t € (or,..., a m ) and 
thus C ai: ...,a m = 0. 

Proof of Lemma S.4: Without any loss of generality, we assume a\ > «2 > 

• • • > • Then maxi<j< m at > J implies a\ > J. In this case, for any 

r € (at,..., a m ), r~ 1 ( 1) := {i \ r(i) = 1} C {1,..., d + J — 1} and |r _1 (l)| = 
ai • If |t -1 (1) n{l,...,d}| > 2, then \A r \ = 0 due to Lemma S.3; otherwise 



d-\-s,d-\-t 


dl dl 

d6 o dOf 
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{d + 1,..., d + J — 1} C r 1 (1) and thus \A r \ = 0 due to Lemma 2.1. Thus 
c Qll ..., Qm = 0 according to (2.3) provided in Theorem 2.2. □ 

Lemma S.5. //#{?' : at > 1} < d, then |Ar| = 0 for any r £ (cg, ..., a m ) and 
thus c aij ... jQ , m = 0. 

Proof of Lemma S.5: Without any loss of generality, we assume a\ > a 2 > 
• • • > a m . Then #{i : oij > 1} < d indicates ot d +\ = • • • = a m = 0. Let 
r : {1, 2,... , d + J — 1} —> {1, —j m} satisfy r £ (cti,..., a m ). Then the (d + 
J — 1) x {d + J — 1) matrix A r can be written as 

( A t i A t 2 \ 

\ A t 3 A t 4 J 

_ j (&T(s)3'T(s)s3'T(s)t)s=l,...d;t=l,...,d ( •£-r(.s).sC'r(s)4).s=l,...,c(;t=l,..., J— 1 

\ ( ^T(d+s)s-^T(d+s)t)s=l,...,J—l;t=l,...,d A T 4 

where the (J — 1) x (J — 1) matrix A T 4 is either a single entry u T ( r j+i)\ (if J = 2) 

or symmetric tri-diagonal with diagonal entries u T ^ d+1 y l ,... ,u T r d+J _iyj_i, upper 

off-diagonal entries — b T i d +\)2, ■ ■ ■, —■ K(d+J-2),j-i, and lower off-diagonal entries 

—b T ( d+ 2 ) 2 ; • • •, —. Note that A r is asymmetric in general. 

If ff{i : oti > 1} < d — 1, then there exists an io such that 1 < io < d and 

|r _1 (i 0 ) H {1,..., d}| > 2. In this case, \A T \ = 0 according to Lemma S.3. 

If ff{i : oti > 1} = d, we may assume |r _1 (z) fl {1,... , d}| = 1 for i = 1,..., d 

(otherwise \A r \ = 0 according to Lemma S.3). Suppose an > a 2 > • ■ • > > 

2 > afc+i . Then {d+1,... ,d+ J— 1} C uf =1 r^ 1 (z) and Yli=i( a i~ 1) = J~ 1- In 

order to show \A T \ = 0, we first replace A r i with A^l = {e T ( s )X T r s ^ t ) s =i,...d-, t=i,...,d 

and replace A r2 with A^ = (—c T ^ t ) s =i t ... :d; t=i,...,J-i- It changes A r into a new 

matrix A^ . Note that \A T \ = fl 1 =1 x t(s)s • lAi- 1 ^. According to Lemma S.2, 

the sum of the columns of A^ is (—e T m,..., —e r (d)) T , and the elementwise sum 

of the columns of A t4 is (c T(d+1)1 , c T(d+2 ) 2 , ■■■, c r(d+ j_ 1))J _ 1 ) T . Secondly, for 

t 1, • • ■, d, we add ®t(i))■***> ®r(d)> Cr(d- 1 - 1)1 > • ■ ■) c r(d+j—l),j-i) to the 

tth column of A^ -. We denote the resulting matrix by Note that |A^| = 
(< 2 \ ( 2 ) 

\A r '|. We consider the sub-matrix A K i which consists of the first d columns of 
A^. For s £ r -1 ( 1), the sth row of A^j is simply 0. For i = 2,..., k, the jth 
row of A^j is proportional to (xn — xn,Xi 2 — x 42 ,... ,Xi d — x\ d ) if j £ r _1 (i). 
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Therefore, Rank(A^) < (d + J — 1) — a\ — Yli=2( a i — 1) = o? — 1, which leads 
to \A^\ = 0 and thus || = 0, \A r \ = 0. According to (2.3) in Theorem 2.2, 

c ai,...,a m d- f3 

Lemma S.6. F = F(p) is always positive semi-definite. It is positive definite if 
and only if p € S+ . Furthermore, log/(p) is concave on S. 

For Section 5.2: The procedure seeking for analytic solutions here follows 
Tong, Volkmer, and Yang (2014). As a direct conclusion of the Karush-Kuhn- 
Tucker conditions (see also Theorem 5.10), a necessary condition for (pi,p 2 -/P:i) to 
maximize f(p\,p 2 ;P:i) in (5.5) is (5.6), which are equivalent to df/dp \ = df /dp 3 , 


df/dp2 = df /dps . In terms of Pi,wfi s, they are 

(P3 - Pi)(pim + P2W2 + P3W3) = (W 3 - U>l)PlP 3 (S.5) 

(P3 ~P2)(p\W\ + P2W2 +P3W3 ) = {W3 - W 2 )P 2 P 3 (S. 6 ) 

Denote y\ = P 1 /P 3 > 0 and y 2 = P 2 /P 3 > 0. Since p± + p 2 + P 3 = 1, it implies 
P3 = l/(?/i + 2/2 + 1), Pi = 2 / 1 /( 2/1 + 2/2 + 1), and p 2 = 2 / 2 /( 2/1 + 2/2 + 1)- In terms 
°f 2 / 1 ; 2 / 2 ; (S.5) and (S.6) are equivalent to 

(1 - 2 / 1 )( 2/1 ■wi + y 2 w 2 + w 3 ) = (w 3 -wi)yi (S.7) 

(1 - 2 / 2 )( 2 / 1 Wl + V2W2 + wfi) = ( w 3 -w 2 )y2 (S.8) 


Lemma S.7. Suppose 0 < w 3 < w 2 < w\ . If (pi,P 2 -P 3 ) maximizes f(pi,P 2 ,P 3 ) 
in (5.5) under the constrains pi,p 2 ,P 3 > 0 and p\ + p 2 + P 3 = 1, then 0 < p 3 < 
P2<Pi < 1- 

The proof of the lemma above is straightforward. Because otherwise one 
could exchange Pi,Pj to strictly improve f(pi,p 2 ,P 3 ) ■ Now we are ready to get 
solutions to equations (S.7) and (S.8) case by case. 

(i) wi = w 3 . In that case, (S.7) implies 2/1 = 1. After plugging it into (S.8), the 
only positive solution is 2/2 = (—3uq + 2 w 2 + y/9wi — 4wiw 2 + 4wf)/ ( 2 w 2 ) ■ 

(ii) w 2 = W 3 . In that case, (S.8) implies 2/2 = 1. After plugging it into (S.7), the 
only positive solution is y\ = ( 2 w\ — 3w 2 + y/— 4wiiv 2 + 9w%)/(2wi) . 
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(iii) w\ = W2 but w\ ^ W3 and W2 / W 3 . The ratio of (S.7) and (S.8) leads 

to yi = y2 ■ After plugging it into (S.7), the only positive solution is 
y\ = (3w\ — 2w3 + a /9 wf — 4w\W3 + /(Aw\) . 

(iv) wi,W 2 ,u >3 are distinct. Without any loss of generality, we assume 0 < W 3 < 
w 2 < W{ . (Otherwise the previous elimination procedure in the order of 
P 3 iP 2 ,Pi could be easily changed accordingly.) Based on Lemma S.7, if 
(pi ■/P 2 ,p:i) maximizes /, 4 , then 0 < P 3 < P 2 < pi < 1 and thus y\ > y 2 > 1. 
The ratio of (S.7) and (S. 8 ) leads to (1 — pi)/(l — 2 / 2 ) = ( 1 C 3 — w\)/(w 3 — 
w 2 ) ■ 2/1 / 2 / 2 , which implies 

(wi - w 3)Vi / Sq n 

2 (w 2 - w 3 ) + (wi - w 2 )yi 

Note that {w2 — W3) + (w\ — W2)yi > w\ — W3 > 0. After plugging (S.9) into 
(S.7), we get 

co + cyyi + c 2 y\ + c 3 pi = 0 (S.10) 

where Co = W 3 (iU 2 — W 3 ) > 0 , c± = “&W\W 2 — W 1 W 3 — A 1 U 2 W 3 + 2 w% > 0, 
C2 = 2 w\ — Aw\W2 — W1W3 + 3W2W3, C3 = wi(w 2 — W\) < 0 . 

Lemma S. 8 . Suppose 0 < W3 < W2 < w± . Then equation (S. 10 ) has one and 
only one solution y\ > 1 . Furthermore, y\> 1 . 


Proof of Lemma S.8: In order to locate the roots of equation (S.10), we let 
fi(yi) = co + cipi + c 2 y\ + c 3 y\. Then /i(l) = c 0 + ci + c 2 + c 3 = (wi-w 3 ) 2 > 0. 

On the other hand, the first derivative of /1 is = ao+«i?/i+« 2 'l/i , where 

Oo = 3wiW 2 — WiW3 — 4w 2 W3 + 2wl = Wi{w2—W3) + 2{wi—W2)W2 + 2{w2—W3) 2 > 0, 
ai = 2(2 w 2 — 4 w\W 2 — W 1 W 3 + 3 W 2 W 3 ), and <22 = 3 wi(vj 2 — vj\) < 0. Therefore, 
af - 4 a 0 a 2 > a{ > 0 and f[(yi) = a 2 (yi - yu)(yi - 2 / 12 ), where 


-ai + y/aj - 4a 0 a 2 ^ _ -ai - y/a\- 4a 0 a 2 ^ 

y n =-V-< vi2 = -V-> v n 


2a2 


2 a 2 


It can be verified that y \ 2 < 1 if and only if w\ < 2 {w 2 + W3). There are two 
cases: Case (i): If P 12 < 1, then f[(yi) < 0 for all y\ > 1. That is, fi(yi) strictly 
decreases after y\ = 1. Since /i(l) > 0 and /i(oo) = —00, then there is one and 
only one solution in (l,oo); Case (ii): If 2/12 > 1, then f[(y\) > 0 for y\ £ [ 1 , 2 / 12 ] 



SUPPLEMENTARY MATERIALS 


S7 


and f[(yi) < 0 for y\ € (7/12,00). That is, fi(yi) increases in [1, 2/12] and then 
strictly decreases in (7/12,00). Again, due to /i(l) > 0 and /i(oo) = —00, there 
is one and only one solution in (1, 00). In either case, the conclusion is justified. 
□ 


S.3 Additional Proofs 


Proof of Theorem 2.1 As a direct conclusion of Lemma S.l and Lemma S.2, 
we obtain Theorem 2.1. □ 

Examples of A *3 in Theorem 2.1 include ( 77 * 1 ), 


f un -b i2 
V - 6*2 Ui 2 


/ un - 6*2 0 

—bi 2 Ui 2 —bis 
\ 0 - 6*3 Ui3 


( Un - 6*2 0 0 ^ 

-6*2 Ui 2 -6*3 0 

0 -6*3 Ui 3 -6*4 

\ 0 0 - 6*4 1t*4 / 


for J = 2, 3, 4, or 5 respectively. 


Proof of Theorem 2.2 To study the structure of |F| as a polynomial function 
of (ni,...,n m ), we denote the (k,l )th entry of A* by a£j. Given a row map 
r : {1,2,..., d+ J — 1} —> {1,..., to}, we define a (d+ J — 1) x (d + J — 1) matrix 
At = (4f») whose fcth row is given by the fcth row of A r ( fc ) . For a power 
index ( 04 ,..., a m ) with a* € {0,1,... , d, + J — 1} and ck* = d + J — 1, we 
denote 

t € (ai, • • •, o m ) 

if a* = #{j : r(j) = 1 } for each i = 1,..., to. In terms of the construction of A r , 
it says that a* rows of A r are from the matrix A* . 

According to the Leibniz formula for the determinant, 


F 


m 

Y n*A* 

i =1 


E (- 1 ) sgn(l 




CT SSd + j_i 


d+J—1 m 

n e rH(j( kMk) 


k =1 i =1 


where cr is a permutation of {1,2,..., d + J — 1}, and sgn(cr) is the sign or 
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signature of a. Therefore, 

d+J—l 

.= E G 1 )' 6 ” 1 '” E II <*$?> 

u&Sd+j -1 re(ai,...,a m ) k=l 

d+J—l 

e e n <$>> 

re(ai,...,a m ) ceSd+j-i k =1 

= E IAI 

Tg(ai,...,a m ) 

□ 


Proof of Lemma 2.2 To simplify the notations, we let i s = s +1, s = 0,..., d. 
That is, ai = J — 1, «2 = • • • = aj+i = 1. There are only two types of 
r € («i,..., a m ), such that, |A r | may not be 0. 

t of type I: There exist 1 < k < d, 2 < l < d + 1, and 1 < q < J — 1, such 
that, t(/c) = 1 and r(d + q) = l. Following a similar procedure as in the proof of 
Lemma S.5, we obtain 

d -\-1 d 

|A r | = e, ■ |A 13 | • (—l) d |Xi[l, 2,..., d + 1]| • (-l) sgn ^ J] x T{s)s • ^ 

i=2 s=l &l 

t of type II: r{d +1) = • • • = r(d + J — 1) = 1 and {r(l),..., r(d)} = 
{2, + 1}. It can be verified that 

d +1 d 

i a t i= n ei • i a 13 i • (-i^ix^i, 2 ,..., d +in • (-i) sgn ^ n x T(s)s 
i =2 s=l 

According to Theorem 2.2, 

c a lt ...,am — y ] + y ] l-A-rl 

t of type I r of type II 

d +1 / d d+1 

= Ilej - |Ais| • (-1)"|X 1 [1,2.rf+ 1]| - EE E 

i=2 \^k- 1 1=2 T€Sd+i‘T(k)=l,T(d-\-l)=l 
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d 

(_l)sgn(r) JJ 

x t ( s)s 

J-l 

E?+ 

E 



3 = 1 


q=l tl 

TeSd+vr(d+l)=l s =1 / 

d+ 1 





d 

ik 

■ |Ai 3 • 

(-l, d 

Xi[l, 2,. 

• • > d + 1] • 


i=2 





reSd +1 a=l 

d+ 1 






ik 

i = 2 

• |Ai 3 • 

(-1)“ 

Xi[l, 2,. 

• • j d + 1] ' 

.(-l) d |X 1 [l,2,...,d+l]| 

d+ 1 






IK 

i=2 

• |Ai 3 • 

|Xi[i 

, 2,..., d + 1] 2 



where S r j + \ is the set of permutations of {l,...,d + l}. The general case with 
*o, *i ,...,id can be obtained similarly. □ 


Proof of Theorem 2.4 Suppose Rank(Xi) = d + 1. Then there exist io> • • •, 
id € {1, such that, |Xi[*o, ii, ■■■, id\ j 7 ^ 0. According to Lemma S.4, 

/(p) can be regarded as an order-( J — 1) polynomial of pi 0 . Let pi 0 = x £ (0,1) 
and pi = (1 — x)/(m — 1) for i ^ iq . Based on Lemma 2.2, /(p) can be written 
as 


fi 0 ( x ) = aj-ix 


J-i 


H-h a\x 


1 — x 
m — 1 
1 — x 
m — 1 


+ aj- 2 x 


J—2 


d+J-2 


+ ao 


1 — X 

m — 1 
1 — x 
m — 1 


d +1 


d+J—1 


where 


aj-i = |A,; o3 | JJe*/ |Xi[i 0 ,i'i,...,^]| 2 > 0 

{ii,...,^}c{ 1 ,...,m}\{i 0 } s=l 


Therefore, lim a ,_ >1 -(l — x)~ d x 1 ~ J fi 0 (x) = (to — l) _d aj_i > 0. That is, /(p) > 0 
for pi 0 = x close enough to 1 and pi = (1 — x)/(m — 1) for i ^ io . 

In order to justify that the condition Rank(Xi) = d + 1 is also necessary, 
we only need to show that /(p) = 0 if Rank(Xi) < d. Actually, for any r : 
{1,..., d+ J — 1} — > {1,..., m}, we construct as in the proof of Lemma S.5. 
Then lAri = nJLi x T ( s ) s ■ jA^I. Similar as in the proof of Lemma S.5, for 
/ 1, . . . , d, we add X T (^t ( j • ■ ■ > ^r(d): Cr(d- (-1)1 > • • • j ^r(d+ J— 1 ), J~ 1 ) fo the 

ith column of A^. We denote the resulting matrix by A®. Note that |A^| = 
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(3) . (3) 

| At |. We consider the sub-matrix which consists of the first d columns of 
A®. For s £ t _1 (t( 1)), the sth row of is simply 0. For s = 2,..., k, the 
sth row of A K J is e T ( s )(x T (s)i - x T (i)i, • • •, x r[s)d - x r (i)d)- For s = 1,..., J - 1, 
the (d + s)th row of A^j is -c T{d+s)s (x T{d+s)1 - ® T (i)i>.. • ,x r{d+s)d - ® T (i)d)- 
We claim that Rank(A^) < d — 1. Otherwise, if Rank(A^) = d, then there 
exist i \,.-.., i d £ {2,..., d + J — 1}, such that, the sub-matrix consisting of the 
iith,..., Oth rows of A^ is nonsingular. Then the sub-matrix consisting of the 
r(l)th, r(?'i)th,..., r(0)th rows of Xi is nonsingular, which implies Rank(Xi) = 

/o\ /ox 

d + 1. The contradiction implies Rank(A^ f/ ) < d — 1. Then | At | = 0 and thus 
|A r | = 0 for each r. Based on Theorem 2.2, |F| = 0 and thus /(p) = 0. □ 

Proof of Theorem 3.5 Combining Theorem 2.1 and Theorem 2.4, it is straight¬ 
forward that /(p) = 0 if Rank( Xi[{i | p % > 0}]) < d. We only need to show 
that /(p) > 0 if Rank(Xi[{i | pi > 0}]) = d + 1. Due to Theorem 2.1, we only 
need to verify the case pt > 0 ,i = 1,... ,m. (Otherwise, we may simply remove 
all support points with p % = 0.) 

Suppose pi > 0 ,i = 1 and Rank(Xi) = d + 1. Then there exist 

i(h • • • iid £ ,?n}, such that, |Xi[*o,... ,*d]| / 0. According to the proof 

of Theorem 2.4, for each i £ {?'o, • • • ,id}, there exists an £ (0,1), such that, 
/(p) > 0 as long as pi = x £ (1 — e*, 1) and pj = (1 — x)/(m — 1) for j / i. On 
the other hand, for each i ^ {io, ..., i d }, if we denote the jth row of Xi by a.j, 
j = 1 ,..., m, then ai = ao«i 0 + • • • + a d oti d for some real numbers ao,...,a d ■ 
Since a t / 0, then at least one a* / 0. Without any loss of generality, we assume 
a o / 0. Then it can be verified that |Xi [i, i \,..., i d \ \ / 0 too. Following the proof 
of Theorem 2.4 again, for such an i £ {io,..., i d }, there also exists an q £ (0,1), 
such that, /(p) > 0 as long as pi = x £ (1 — e*, 1) and pj = (1 — x)/(m — 1) for 
j 7 ^ i. Let e* = niin{mirije^, (m — l)min,pi, 1 —l/m}/2. For i = 1,...,m, denote 
= (5?;i, • • •, 5i m ) T £ S with §a = 1 — e* and Sij = e */(m — 1) for j ^ i. It can be 
verified that p = ai<5i + - • ■ + a m 8 m with a* = (pi — e*/(m —1))/(1 — 1)). 

By the choice of e*, f(5i) > 0, > 0, i = 1,..., rn. and Yli a i = 1- Then /(p) > 0 

according to Lemma S.6. □ 

Proof of Corollary 5.3 In order to check when a minimally supported design 
supported only on {xi,X 2 } is D-optimal, we add one more support point, that 
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is, x 3 . According to Theorem 2.2, Lemmas S.4, S.5, and 2.2, the objective 
function for a D-optimal approximate design on {at, x 2 , £ 3 } is f(pi,P2,P3) = 
PlP2 (c210.Pl + C120P2) +Pl/>3(C201P1 + C102P3) + P 2 P 3 (C 021 P 2 + C012P3) + CmPlP 2 P 3 , 
where 

C 210 = e 2 5ii5i2( 7r ii 7r i27ri 3 )' 1 (a;i - x 2 ) 2 > 0 

C120 = ei£f| 1 £f| 2 (7r2l7r 2 27r23)' 1 (xi - x 2 ) 2 > 0 

C 201 = e 3 gj l g‘f 2 ('Kn^i2^i3)~ 1 (xi - x 3 ) 2 > 0 

C102 = ei5| 1 £f| 2 ( 7r 3l7!'327r33) _1 (®l — X 3 ) 2 > 0 

C021 = e35| 1 £f| 2 (7r2l7r227T2 3 )' 1 (x2 - A 3 ) 2 > 0 

C012 = e 2 fifil 5 ' 32 ( 7 r 3 l 7 r 32 Vr 3 3 )' 1 (x 2 ~x 3 ) 2 > 0 

Cm = e 1 (u 22 u 31 + u 21 u 32 - 2 b 22 b 3 2 )(x 1 - x 2 )(x 1 - x 3 ) + 

e 2 (u 12 u 31 + U11U32 - 2 b V 2 b 32 )(x 2 - x 1 )(x 2 - x 3 ) + 

e 3 (u 12 u 2 i + u u u 22 - 2 b 12 b 22 )(x 3 - x{)(x 3 - x 2 ) 

Based on Theorem 5.10, the design p = (PiPP 2 ,0) T is D-optimal if and only if 

df{p)/df(p 1 ) = df{p)/df(p 2 ) > df{p)/df(p 3 ) 

Similar conclusions could be justified for X 4 ,..., x rn if to > 4. □ 


Proof of Theorem 5.12 According to the solutions provided by the software 
Mathematica, the largest root of equation (S.10) after simplification is 


b 2 2 1 / 3 (36 1 - b%) A 1 / 3 

y 3AV3 + 3 x 2V3 


(S.ll) 


where A = -27 b 0 + 96i& 2 - 26 3 + 3 3 / 2 (276 2 + 46 3 - 186 0 6i6 2 - b\b\ + 46 0 fe|) 1/2 , 
and bi = Ci/c 3 , i = 0,1, 2. Note that the calculation of A and thus y\ should be 
regarded as operations among complex numbers since the expression under square 
root could be negative. Nevertheless, y\ at the end would be a real number. Thus 
we are able to provide the analytic solution maximizing f{pi,p 2 ,p 3 )- □ 


Proof of Corollary 5.5 In order to check when a minimally supported design 
is D-optimal, we first add the four design points, that is, we consider four design 



S12 


JIE YANG, LIPING TONG AND ABHYUDAY MANDAL 


points (xn,Xi 2 ), i = 1,2, 3,4 and check when the D-optimal design could be 
constructed on the first three design points. Let Xi be defined as in Lemma 2.2. 
In this case, Xi is a 4 x 3 matrix. Following Theorem 2.2, Lemmas S.4, S.5, and 
2.2, the objective function for a minimally supported design at (d, J, m) = (2, 3,4) 
is 


f(Pl,P2,P3,P4) 


= CllUPlP2P3P4 

+ |Xi[l, 2,3]| 2 eie 2 e 3 • piP2P3(wiPi + w 2 p 2 + w 3 p 3 ) 
+ |Xi[l,2,4]| 2 eie 2 e 4 • PiP2 P4 (wiPi + w 2 p - 2 + mp4) 
+ |Xi[l,3,4]| 2 eie 3 e 4 • PiP3P4(wiPi + w 3 p 3 + w 4 p 4 ) 
+ |Xi[2,3,4]| 2 e 2 e 3 e± ■ p 2 p 3 Pi{v)2P2 + w 3 p 3 + w 4 p 4 ) 


where e % = u lA + u i2 - 2 b i2 , = e i 1 1 , i = 1 , 2 ,3,4, and 


C 1111 = ^2 e i e Mki u i2 + u k2 un - 2b k2 bi 2 ) ■ |X a [i, j. A:] | • |X 4 [?:, j, Z]| (S.12) 

l<i<j<4 

with {i, j. k, 1} = {1, 2,3,4} given 1 < i < j < 4. 

According to Theorem 5.10, a minimally supported design p = (pl,p 2 ,P 3 , 
0) T in this case is D-optimal if and only if df /dpi = df/dp 2 = df/dp 3 > df /dpi 
at p. Then df /dpi = df /dp 2 = df/dp 3 at p is equivalent to (1) of Corollary 5.5, 
and df /dpi < df /dpi at p leads to (2) of Corollary 5.5 since the forms of df/dpi 
at p, i = 1, 2, 3 will not change if more than four design points (i.e., m > 4) are 
added into consideration. Note that 

|Xi[l, 2, 2>\\ 2 eie 2 e 3 p* 2 p 3 (2wip\ + w 2 p* 2 + w 3 p^) in (2) of Corollary 5.5 is equal to 
df/dpi at p. It could be replaced with |Xi[l, 2, 3\\ 2 eie 2 e 3 p\p% i (wip\ + 2w 2 p 2 + 
w 3 p* 3 ) (i.e., df/dp 2 ), or |Xi[l, 2, ^eie 2 e 3 p\pl{wip\+w 2 p* 2 +2w 3 pf) (i.e., df/dp 3 ), 
since these three are all equal. □ 


S.4 Maximization of fi(z) in Section 3 

According to Theorem 3.6, fi(z) is an order- [d + J — 1) polynomial of z. In other 
to determine its coefficients ao, ai ,..., aj_i as in (3.2), we need to calculate /*(0), 
/i(l/2), fi( 1/3), ..., fi(l/J), which are J determinants defined in (3.1). 

Note that Bj^ is a matrix determined by J — 1 only. For example, Bf 1 = 1 
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for J = 2, 



for J = 3, 4, or 5 respectively. 

Once ao,...,aj-i in (3.2) are determined, the maximization of fi(z) on 
z G [ 0 , 1 ] is numerically straightforward since it is a polynomial and its derivative 
is given by 

J-i J-i 

3 =1 3=0 

(S.13) 

S.5 Exchange algorithm for D-optimal exact allocation in Sec¬ 
tion 4 

Exchange algorithm for D-optimal allocation (n ±,..., n m ) T given n > 0: 

1° Start with an initial design n = (m,..., n m ) T such that /(n) > 0. 

2° Set up a random order of (i.j) going through all pairs {(1,2), (1,3), ..., 
(1, m), (2,3), ■ ■ •, (m - l,m)}. 

3° For each let c = ni + nj . If c = 0, let n*- = n. Otherwise, there are 

two cases. Case one: 0 < c < J, we calculate fij(z) as defined in (4.1) for 
z = 0,1,..., c directly and find z* which maximizes fij(z). Case two: c> J, 
we first calculate fij(z) for z = 0,1,..., J; secondly determine Co, ci,..., cj 
in (4.2) according to Theorem 4.9; thirdly calculate fij(z) for z = J+l,..., c 
based on (4.2); fourthly find z* maximizing f tJ (z) for z = 0,... ,c. For both 
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cases, we define 

n ij = (h\, • ■ • j Hi— 1 , Z , Tli- fl, • • • , Hj— 1 , C Z , • • • , ) 

Note that /( n*-) = fij(z*) > /( n) > 0. If /( nT) > /(n), replace n with 
n *j , and /(n) with /( n?-). 

4° Repeat 2° ~ 3° until convergence, that is, /( n*) = /( n) in step 3° for any 

ihj)- 

S.6 More on Example 4.6: Polysilicon Deposition Study 

Table S.l shows the list of experimental settings for the polysilicon deposi¬ 
tion study. The factors are decomposition temperature (d), decomposition 
pressure(R), nitrogen flow (C), silane flow(H), setting time(E), cleaning 
method(I ? ). Column 1 provides original indices of experimental settings out of 
729 distinct ones. For each experimental setting labelled “1” in a design, 9 re¬ 
sponses are collected (Phadke, 1989) and assumed to be independent. 

Table S.2 shows the model matrix for the D-optimal design n Q found for the 
polysilicon deposition study. In this table, each 3-level factor is represented by its 
linear component and quadratic component. Thus there are level combinations 
of 12 predictors. 
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Table S.l: Polysilicon Deposition Study: Experimental Settings for the Original Design, 
Rounded Approximate Design and D-optimal Exact Design 


Index 

A 

B 

C 

D 

E 

F 

Original 

Rounded 

D-optimal 

1 

1 

1 

1 

1 

1 

1 

1 

0 

0 

76 

1 

1 

3 

3 

2 

1 

1 

0 

0 

89 

1 

2 

1 

1 

3 

2 

1 

0 

0 

98 

1 

2 

1 

2 

3 

2 

0 

0 

1 

111 

1 

2 

2 

1 

1 

3 

0 

0 

1 

116 

1 

2 

2 

1 

3 

2 

0 

1 

0 

122 

1 

2 

2 

2 

2 

2 

1 

0 

0 

130 

1 

2 

2 

3 

2 

1 

0 

0 

1 

167 

1 

3 

1 

1 

2 

2 

0 

0 

1 

181 

1 

3 

1 

3 

1 

1 

0 

1 

0 

199 

1 

3 

2 

2 

1 

1 

0 

1 

1 

201 

1 

3 

2 

2 

1 

3 

1 

0 

0 

243 

1 

3 

3 

3 

3 

3 

1 

0 

1 

258 

2 

1 

1 

2 

2 

3 

1 

0 

0 

286 

2 

1 

2 

2 

3 

1 

0 

1 

0 

290 

2 

1 

2 

3 

1 

2 

1 

0 

0 

291 

2 

1 

2 

3 

1 

3 

0 

1 

0 

294 

2 

1 

2 

3 

2 

3 

0 

0 

1 

299 

2 

1 

3 

1 

1 

2 

0 

0 

1 

301 

2 

1 

3 

1 

2 

1 

0 

1 

0 

313 

2 

1 

3 

2 

3 

1 

0 

0 

1 

331 

2 

2 

1 

1 

3 

1 

0 

1 

1 

336 

2 

2 

1 

2 

1 

3 

0 

1 

1 

339 

2 

2 

1 

2 

2 

3 

0 

1 

0 

350 

2 

2 

1 

3 

3 

2 

0 

1 

0 

365 

2 

2 

2 

2 

2 

2 

0 

0 

1 

376 

2 

2 
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